/************************/
/**    EQUILIBRIUM    ***/
/************************/

double ComputeEquilibrium(double *PARAM){

    int iter, iterprice, igrid, tgrid, ygrid, igridCOH, iteration, iteration2;
    double axopt, *genmoments, GMMmin;
    genmoments = (double *) calloc((nbmoments), sizeof(double));

    // STATISTICS
    double capitalin, capitalout, laborin, laborout, taxoutSEA, taxin, taxout, critprice, critparams, taxrateout;
    
#if OPT_GE == 0
    rstart   = 0.016510391666088;
    lstart   = 1.094347964821837;
    taxstart = 0.009027144467390;
    wstart   = 0.253706723228092;
    capstart = 3.850972708693789;
    #if Teq == 1
        taxstart = 0.009314127952320;
    #endif
    #if Weq == 1
        wstart = 0.253781277852715;
    #endif
    #if Req == 1
        rstart = 0.016489214609214;
    #endif
#endif
    
    /******* ENDOGENOUS PARAMETERS *******/
    betapar     = PARAM[0];             // discount factor
    kappaE      = PARAM[1];             // matching param
    kappaW      = PARAM[2];             // matching param function worker
    psiparW     = PARAM[3];             // min val
    psiparE     = psiparW;
    nupar       = PARAM[4];             // borrowing constraint
    court_cost  = PARAM[5];             // return to scale param
    g1          = PARAM[6];             // mapping 1
    g2          = PARAM[7];             // mapping 2
    g3          = PARAM[8];             // mapping 3
    sigz        = PARAM[9];             // variance of z shock
    pz          = PARAM[10];            // persistence of z shock
    gmapping[0] = g1;
    gmapping[1] = g2;
    gmapping[2] = g3;

        
    // TRANSITION + STATE Z OF ENTREPRENEURS //
    tauchenfun(pz, 1.0, 1.0, sigz, statez, mzprob);
    inv_distri(mzinv, mzprob);
    double meanz = 0.0;
    for(int zgrid = 0; zgrid < maxfirmtype; zgrid++){
        //statez[zgrid] = exp(statez[zgrid]);
        meanz += statez[zgrid]*mzinv[zgrid];
    }
    for(int zgrid = 0; zgrid < maxfirmtype; zgrid++){
        statez[zgrid] = (statez[zgrid])/meanz;
    }
    
    // productivite y
    if(maxyprod > 1){
        printf("Transitory Y process: \t");
        tauchenfun(py, 1.0, 0.0, sigy, yprod, myprod);
        inv_distri(invmyprod, myprod);
        for(ygrid=0;ygrid<maxyprod;ygrid++){yprod[ygrid] = exp(yprod[ygrid]);printf("%f \t",yprod[ygrid]);} printf("\n");
        for(ygrid=0;ygrid<maxyprod;ygrid++){printf("%f \t",invmyprod[ygrid]);} printf("\n");
        for(ygrid=0;ygrid<maxyprod;ygrid++){
            for(int ygridn=0;ygridn<maxyprod;ygridn++){
                printf("%f \t",myprod[ygrid][ygridn]);
            } printf("\n");
        } printf("\n");
    }else{
        yprod[0] = 1.0; invmyprod[0] = 1.0; myprod[0][0] = 1.0;
    }

    
    
    /*******************************/
    /******* GUESS ON PRICES *******/
    /*******************************/
    // Update prices if allowed.
    double Rannual, RquartG, Rcomparatif, Wcomparatif;
    // Annual gross interest rate
    // --> quarter gross interest rate
    // --> quarter net interest rate

    rstar       = rstart;
    lstar       = lstart;
    taxrate     = taxstart;
    capitalin   = capstart;
    wstar       = wstart;

        
    
    // GENERATE THE GRID FOR COH //
        
    // compute the optimal capital level -- this define the maximum possible grid level //+- TO CHECK
//    #if OPT_labor_demand == 0
//        axopt  = min(pow((rstar + wedgerate + deltapar)/(statez[(maxfirmtype-1)]*nupar*gmapping[(maxtheta-1)]), (1/(nupar-1))), (1+lambdapar)*(Gridmax));
//    #endif
//        
//    #if OPT_labor_demand == 1
//        axopt  = max(pow((rstar + wedgerate + deltapar)/(statez[(maxfirmtype-1)]*nupar*gammapar*gmapping[(maxtheta-1)]), (1/(nupar*gammapar-1))),pow(((rstar + wedgerate + deltapar)/(statez[(maxfirmtype-1)]*gmapping[(maxtheta-1)]*pow((wstar)/(nupar*statez[(maxfirmtype-1)]*gmapping[(maxtheta-1)]*(1.0-gammapar)),((1.0-gammapar)*nupar)/((1.0-gammapar)*nupar-1.0))))*((1.0-(1.0-gammapar)*nupar)/(gammapar*nupar)), (1/((gammapar*nupar)/(1.0-(1.0-gammapar)*nupar) - 1.0))));
//        axopt  = min((1+lambdapar)*(Gridmax),axopt);
//        printf("\n \n %f %f %f %f \n \n",pow((rstar + wedgerate + deltapar)/(statez[(maxfirmtype-1)]*nupar*gammapar*gmapping[(maxtheta-1)]), (1/(nupar*gammapar-1))),(1+lambdapar)*(Gridmax),pow(((rstar + wedgerate + deltapar)/(statez[(maxfirmtype-1)]*gmapping[(maxtheta-1)]*pow((wstar)/(nupar*statez[(maxfirmtype-1)]*gmapping[(maxtheta-1)]*(1.0-gammapar)),((1.0-gammapar)*nupar)/((1.0-gammapar)*nupar-1.0))))*((1.0-(1.0-gammapar)*nupar)/(gammapar*nupar)), (1/((gammapar*nupar)/(1.0-(1.0-gammapar)*nupar) - 1.0))),axopt);
//    #endif
        
        
    
//    GridmaxCOH = profitE(statez[(maxfirmtype-1)],gmapping[(maxtheta-1)],axopt) + Gridmax + (rstar+wedgerate)*(Gridmax - axopt);
    GridmaxCOH = Gridmax + 200;
//    #if SEA == 1
//        GridmaxCOH = GridmaxCOH + bstarfun(profitE(statez[(maxfirmtype-1)],gmapping[(maxtheta-1)],axopt),(maxtheta-1));
//    #endif

    // grid of COH
    for(igridCOH=0; igridCOH < maxgridCOH; igridCOH++){gridCOH[igridCOH] = expspace(igridCOH,GridminCOH,GridmaxCOH,EchelleCOH,maxgridCOH);}
    
    
    
    #if SEA == 1
        lumpsumSEA = lumpsumSEAstart;   // if financed by labor income tax, does not care of it.
    #endif


        
    /****************/
    /** TEMP FILE  **/
    /****************/

    time_t rawtime,timeofstart,timeofend;
    struct tm * timeinfo;

    char buffer [80];

    FILE    *tempfileoutfile;

    //stop watch
    time ( &rawtime );
    time (&timeofstart);
    timeinfo = localtime ( &rawtime );
    snprintf(buffer, sizeof(char) * 80, "OUTPUT/temp/output_policy=%d_tax=%d_%d_%d_labor=%d_f=%f.out",OPT_policy_exp,Teq,Weq,Req,OPT_labor_demand,DRIpar);
    strncpy(tempfileout, buffer, sizeof(tempfileout) - 1);
    tempfileout[sizeof(tempfileout) - 1] = 0;


    tempfileoutfile=fopen(tempfileout, "w");
    setbuf ( tempfileoutfile , NULL );
    fprintf(tempfileoutfile,"===================Program Starting: %s\n",asctime(timeinfo));
    fprintf(tempfileoutfile,"------ Estimated Model parameters------\n");
    fprintf(tempfileoutfile,"%20s\t\t%20.15f\n","betapar",betapar);
    fprintf(tempfileoutfile,"%20s\t\t%20.15f\t%20.15f\t%20.15f\n","Etapar",etapar[0],etapar[1],etapar[2]);
    fprintf(tempfileoutfile,"%20s\t\t%20.15f\n","Kappa",kappaW);
    fprintf(tempfileoutfile,"%20s\t\t%20.15f\t%20.15f\t%20.15f\n","PROD",prod[0], prod[1],prod[2]);
    fprintf(tempfileoutfile,"%20s\t\t%20.15f\t%20.15f\t%20.15f\n","gmapping",gmapping[0], gmapping[1], gmapping[2]);
    fprintf(tempfileoutfile,"%20s\t\t%20.15f\t%20.15f\t%20.15f\n","gmapping",gmapping[0]/meang, gmapping[1]/meang, gmapping[2]/meang);
    fprintf(tempfileoutfile,"%20s\t\t%20.15f\n","meang",meang);
    fprintf(tempfileoutfile,"%20s\t\t%20.15f\n","kappaE",kappaE);
    fprintf(tempfileoutfile,"%20s\t\t%20.15f\n","psiparW",psiparW);
    fprintf(tempfileoutfile,"%20s\t\t%20.15f\n","psiparE",psiparE);
    fprintf(tempfileoutfile,"%20s\t\t%20.15f\n","nupar",nupar);
    fprintf(tempfileoutfile,"%20s\t\t%20.15f\n","lambdapar",lambdapar);
    fprintf(tempfileoutfile,"%20s\t\t%20.15f\n","phipar",phipar);
    fprintf(tempfileoutfile,"%20s\t\t%20.15f\n","minpar",minpar);
    fprintf(tempfileoutfile,"%20s\t\t%20.15f\n","DRIpar",DRIpar);
    fprintf(tempfileoutfile,"%20s\t\t%20.15f\n","mupar",mupar);
    fprintf(tempfileoutfile,"%20s\t\t%20.15f\n","rhostar",rhostar);
    fprintf(tempfileoutfile,"%20s\t\t%20.15f\n","pLTpar",pLTpar);
    fprintf(tempfileoutfile,"%20s\t\t%20.15f\n","muz",muz);
    fprintf(tempfileoutfile,"%20s\t\t%20.15f\n","pz",pz);
    fprintf(tempfileoutfile,"%20s\t\t%20.15f\n","sigz",sigz);

    // TRANSITION MATRIX Z //
    for(int zgrid = 0; zgrid < maxfirmtype; zgrid++) {
    fprintf(tempfileoutfile,"%20s\t%5d\t","PiZ",zgrid);
        for(int zgrid2 = 0; zgrid2 < maxfirmtype; zgrid2++) {
            fprintf(tempfileoutfile, "%20.15f\t", mzprob[zgrid][zgrid2]);
        }
    fprintf(tempfileoutfile, "\n");
    }
    // STATE Z //
    fprintf(tempfileoutfile,"%20s\t","Zgrid");
    for(int zgrid = 0; zgrid < maxfirmtype; zgrid++) {
        fprintf(tempfileoutfile,"%20.15f\t",statez[zgrid]);
    }
    fprintf(tempfileoutfile, "\n");
        
    // TRANSITION MATRIX Z //
    for(int zgrid = 0; zgrid < maxfirmtype; zgrid++) {
    printf("%20s\t%5d\t","PiZ",zgrid);
        for(int zgrid2 = 0; zgrid2 < maxfirmtype; zgrid2++) {
            printf("%20.15f\t", mzprob[zgrid][zgrid2]);
        }
    printf("\n");
    }
    // STATE Z //
    printf("%20s\t","Zgrid");
    for(int zgrid = 0; zgrid < maxfirmtype; zgrid++) {
        printf("%20.15f\t",statez[zgrid]);
    }
    printf("\n");

    #if SEA == 1
    fprintf(tempfileoutfile,"%20s\t\t%20.15f\n","mupar",mupar);
    fprintf(tempfileoutfile,"%20s\t\t%20.15f\n","DRIpar",DRIpar);
    #endif
    fprintf(tempfileoutfile,"\n");
    fclose(tempfileoutfile);



    /*************************************/
    // INITIALIZE

    double *distUiE, *distUnE, *distUiNE, *distUnNE, *distWWNE, *distWWE, *distJE, *distJNE;
    // worker
    distWWE     = (double *) calloc((ifulldimW), sizeof(double));
    distWWNE    = (double *) calloc((ifulldimW), sizeof(double));
    // unemployed
    distUiE     = (double *) calloc((ifulldim), sizeof(double));
    distUnE     = (double *) calloc((ifulldim), sizeof(double));
    distUiNE    = (double *) calloc((ifulldim), sizeof(double));
    distUnNE    = (double *) calloc((ifulldim), sizeof(double));
    // entrepreneur
    distJNE     = (double *) calloc((ifulldimE), sizeof(double));
    distJE      = (double *) calloc((ifulldimE), sizeof(double));
    #if SEA == 1
        double *distJEi, *distJNEi;
        distJNEi    = (double *) calloc((ifulldimE), sizeof(double));
        distJEi     = (double *) calloc((ifulldimE), sizeof(double));
    #endif


    // ENTREPRENEUR
    double *valueJNE, *collateralJNE, *sloanJNE, *rateJNE, *saveJNE, *C_JNE, *searchJNE_W, *defaultJNE;  //  non excluded  after knowning type zeta
    double *valueJE, *collateralJE,*saveJE,*C_JE, *searchJE_W;  // excluded after knowning type zeta
    valueJNE        = (double *) calloc((ifulldimE), sizeof(double));
    valueJE         = (double *) calloc((ifulldimE), sizeof(double));
    collateralJNE   = (double *) calloc((ifulldimE), sizeof(double));
    collateralJE    = (double *) calloc((ifulldimE), sizeof(double));
    sloanJNE        = (double *) calloc((ifulldimE), sizeof(double));
    rateJNE         = (double *) calloc((ifulldimE), sizeof(double));
    defaultJNE      = (double *) calloc((ifulldimEE), sizeof(double));
    saveJNE         = (double *) calloc((ifulldimEE), sizeof(double));
    saveJE          = (double *) calloc((ifulldimEE), sizeof(double));
    C_JNE         = (double *) calloc((ifulldimEE), sizeof(double));
    C_JE          = (double *) calloc((ifulldimEE), sizeof(double));
    searchJE_W      = (double *) calloc((ifulldimEE), sizeof(double));
    searchJNE_W     = (double *) calloc((ifulldimEE), sizeof(double));
    #if SEA == 1
    double *valueJNEi, *collateralJNEi, *sloanJNEi, *rateJNEi, *saveJNEi, *C_JNEi, *searchJNEi_W, *defaultJNEi;  // INSURED non excluded after knowning type zeta
    double *valueJEi, *collateralJEi,  *saveJEi, *C_JEi, *searchJEi_W;  // INSURED excluded  after knowning type zeta
    valueJNEi        = (double *) calloc((ifulldimE), sizeof(double));
    valueJEi         = (double *) calloc((ifulldimE), sizeof(double));
    collateralJNEi   = (double *) calloc((ifulldimE), sizeof(double));
    collateralJEi    = (double *) calloc((ifulldimE), sizeof(double));
    sloanJNEi        = (double *) calloc((ifulldimE), sizeof(double));
    rateJNEi         = (double *) calloc((ifulldimE), sizeof(double));
    defaultJNEi      = (double *) calloc((ifulldimEE), sizeof(double));
    saveJNEi         = (double *) calloc((ifulldimEE), sizeof(double));
    saveJEi          = (double *) calloc((ifulldimEE), sizeof(double));
    C_JNEi           = (double *) calloc((ifulldimEE), sizeof(double));
    C_JEi            = (double *) calloc((ifulldimEE), sizeof(double));
    searchJEi_W      = (double *) calloc((ifulldimEE), sizeof(double));
    searchJNEi_W     = (double *) calloc((ifulldimEE), sizeof(double));
    #endif

        
    // WORKER
    double *valueWWE, *saveWWE, *C_WWE,  *searchWWE_J, *valueWWNE, *saveWWNE, *C_WWNE, *searchWWNE_J;
    valueWWE        = (double *) calloc((ifulldimW), sizeof(double));
    valueWWNE       = (double *) calloc((ifulldimW), sizeof(double));
    saveWWE         = (double *) calloc((ifulldimW), sizeof(double));
    saveWWNE        = (double *) calloc((ifulldimW), sizeof(double));
    C_WWE           = (double *) calloc((ifulldimW), sizeof(double));
    C_WWNE          = (double *) calloc((ifulldimW), sizeof(double));
    searchWWE_J     = (double *) calloc((ifulldimW), sizeof(double));
    searchWWNE_J    = (double *) calloc((ifulldimW), sizeof(double));

        
    // UNEMPLOYED
    double *valueUiE, *saveUiE, *C_UiE, *searchUiE_J, *searchUiE_W, *valueUiNE, *saveUiNE, *C_UiNE, *searchUiNE_J, *searchUiNE_W;    /* insured */
    double *valueUnE, *saveUnE, *C_UnE, *searchUnE_J, *searchUnE_W, *valueUnNE, *saveUnNE, *C_UnNE, *searchUnNE_J, *searchUnNE_W;    /* not insured */
    valueUiE     = (double *) calloc((ifulldim), sizeof(double));
    valueUnE     = (double *) calloc((ifulldim), sizeof(double));
    valueUiNE    = (double *) calloc((ifulldim), sizeof(double));
    valueUnNE    = (double *) calloc((ifulldim), sizeof(double));
    saveUiE      = (double *) calloc((ifulldim), sizeof(double));
    saveUnE      = (double *) calloc((ifulldim), sizeof(double));
    saveUiNE     = (double *) calloc((ifulldim), sizeof(double));
    saveUnNE     = (double *) calloc((ifulldim), sizeof(double));
    C_UiE        = (double *) calloc((ifulldim), sizeof(double));
    C_UnE        = (double *) calloc((ifulldim), sizeof(double));
    C_UiNE       = (double *) calloc((ifulldim), sizeof(double));
    C_UnNE       = (double *) calloc((ifulldim), sizeof(double));
    searchUiE_W  = (double *) calloc((ifulldim), sizeof(double));
    searchUnE_W  = (double *) calloc((ifulldim), sizeof(double));
    searchUiE_J  = (double *) calloc((ifulldim), sizeof(double));
    searchUnE_J  = (double *) calloc((ifulldim), sizeof(double));
    searchUiNE_W = (double *) calloc((ifulldim), sizeof(double));
    searchUnNE_W = (double *) calloc((ifulldim), sizeof(double));
    searchUiNE_J = (double *) calloc((ifulldim), sizeof(double));
    searchUnNE_J = (double *) calloc((ifulldim), sizeof(double));


        
    //* --- INITIALIZE THE DISTIRBUTION AND VALUE
        
    // initial condition
    distWWNE[inxW(0, 0, 3)] = s_distWWNE;


    // GENERATE VALUES & GRID //
    for(igrid=0;igrid<maxgrid;igrid++){
        grid[igrid] = expspace(igrid,Gridmin,Gridmax,Echelle1,maxgrid);
        
        for(int tgrid = 0.0; tgrid < maxtheta; tgrid ++){
            for(int zgrid = 0.0; zgrid < maxfirmtype; zgrid++) {
            
                valueJE[inxE(igrid,tgrid,zgrid)]    =   0;
                valueJNE[inxE(igrid,tgrid,zgrid)]   =   0;
                
                #if SEA == 1
                valueJEi[inxE(igrid,tgrid,zgrid)]   =   0;
                valueJNEi[inxE(igrid,tgrid,zgrid)]  =   0;
                #endif
            }
            
            for(int ygrid = 0.0; ygrid < maxyprod; ygrid++) {
                valueWWE[inxW(igrid,tgrid, ygrid)]  =   0;
                valueWWNE[inxW(igrid,tgrid, ygrid)] =   0;
            }
            
            valueUiE[inx(igrid,tgrid)]  =   0;
            valueUiNE[inx(igrid,tgrid)] =   0;
            valueUnE[inx(igrid,tgrid)]  =   0;
            valueUnNE[inx(igrid,tgrid)] =   0;
        }
    }




    /** OPTIMAL SEARCHING EFFORT **/
    compute_search();



    /*** PROJECTED METHOD ***/
    #if indexPM == 1
        PM(nbTperiodPM);
        printf("PROJECTED METHOD FINISHED"); getchar();
    #endif
        


    /*** STARTING LOOP OVER PRICES **/

    critprice = 1.0; iterprice = 0.0;
    

    
    // guesses
    //RquartG = alphapar*TFPpar*pow((lstar/capstart),(1.0-alphapar));
    rstar   = rstart; //RquartG - deltapar - wedgerate;             // rate r^d
    wstar   = wstart; //(1.0-alphapar)*TFPpar*pow((capstart/lstar),(alphapar));
    

#if OPT_GE == 1
    while(critprice > epsilonprice){
#endif
        printf("\n ITER ON PRICES:: R=%f, W=%f, K=%f, L=%f, T=%f || Critere: %f \n", rstar, wstar, capitalin, lstar, taxrate, critprice);
        
        // If bug occurs //
        // if (betapar*(1.0+rstar)>=1.1 || rstar < 0.0005 || capitalin < 0.0 || iterprice > itermax)
        if (rstar < 0.0005 || capitalin < 0.0 || iterprice > itermax){
            
            printf("beta condition %20.15f\t%20.15f\n",rstar,betapar*(1.0+rstar));
            critprice = 0.0;
            for(int l = 0.0; l < nbmoments; l++){genmoments[l] = -1000;}
            
        } else {

            /******* VFI Computation *******/
            // TIME //
            time_t rawtime,timeofstart,timeofend;
            struct tm * timeinfo;
            time ( &rawtime );
            time (&timeofstart);
            timeinfo = localtime ( &rawtime );

            #if indexPM == 0
                #if NOPOLICY == 1
                    VFI(valueUiNE, saveUiNE, C_UiNE,   searchUiNE_W, searchUiNE_J,
                        valueUiE,  saveUiE,  C_UiE,    searchUiE_W,  searchUiE_J,
                        valueUnNE, saveUnNE, C_UnNE,   searchUnNE_W, searchUnNE_J,
                        valueUnE,  saveUnE,  C_UnE,    searchUnE_W,  searchUnE_J,
                        valueWWNE, saveWWNE, C_WWNE,   searchWWNE_J,
                        valueWWE,  saveWWE,  C_WWE,    searchWWE_J,
                        valueJE,   saveJE,   C_JE,     searchJE_W,  collateralJE,
                        valueJNE,  saveJNE,  C_JNE,    searchJNE_W, collateralJNE, sloanJNE, rateJNE, defaultJNE);
                #endif
                #if SEA == 1
                    VFI(valueUiNE, saveUiNE, C_UiNE,   searchUiNE_W, searchUiNE_J,
                        valueUiE,  saveUiE,  C_UiE,    searchUiE_W,  searchUiE_J,
                        valueUnNE, saveUnNE, C_UnNE,   searchUnNE_W, searchUnNE_J,
                        valueUnE,  saveUnE,  C_UnE,    searchUnE_W,  searchUnE_J,
                        valueWWNE, saveWWNE, C_WWNE,   searchWWNE_J,
                        valueWWE,  saveWWE,  C_WWE,    searchWWE_J,
                        valueJE,   saveJE,   C_JE,     searchJE_W,  collateralJE,
                        valueJNE,  saveJNE,  C_JNE,    searchJNE_W, collateralJNE, sloanJNE, rateJNE, defaultJNE,
                        valueJEi,  saveJEi,  C_JEi,    searchJEi_W, collateralJEi,
                        valueJNEi, saveJNEi, C_JNEi,   searchJNEi_W,collateralJNEi, sloanJNEi, rateJNEi, defaultJNEi);
                #endif
            #endif
            
            //getchar();
            
            time ( &rawtime );
            time (&timeofend);
            timeinfo = localtime ( &rawtime );


            printf("=================== VFI ended: %s in %f seconds ================= \n",asctime(timeinfo),difftime(timeofend,timeofstart));


            
            timeinfo = localtime ( &rawtime );
            
            /*******  SIMULATION PART  ********/
            #if indexPM == 0
                #if NOPOLICY == 1
                SIMULATION(valueUiNE, C_UiNE, saveUiNE,searchUiNE_W, searchUiNE_J,
                           valueUiE,  C_UiE, saveUiE, searchUiE_W,  searchUiE_J,
                           valueUnNE, C_UnNE, saveUnNE,searchUnNE_W, searchUnNE_J,
                           valueUnE,  C_UnE, saveUnE, searchUnE_W,  searchUnE_J,
                           valueWWNE, C_WWNE, saveWWNE,searchWWNE_J,
                           valueWWE,  C_WWE, saveWWE, searchWWE_J,
                           valueJE,   C_JE, saveJE,  searchJE_W,  collateralJE,
                           valueJNE,  C_JNE, saveJNE, searchJNE_W, collateralJNE, sloanJNE, rateJNE, defaultJNE,
                           &capitalout, &laborout, &taxout, genmoments,
                           distJE, distJNE, distWWNE, distWWE, distUnE, distUnNE, distUiNE, distUiE);
                #endif
                #if SEA == 1
                SIMULATION(valueUiNE, C_UiNE, saveUiNE, searchUiNE_W,  searchUiNE_J,
                           valueUiE,  C_UiE, saveUiE,  searchUiE_W,   searchUiE_J,
                           valueUnNE, C_UnNE, saveUnNE, searchUnNE_W,  searchUnNE_J,
                           valueUnE,  C_UnE, saveUnE,  searchUnE_W,   searchUnE_J,
                           valueWWNE, C_WWNE, saveWWNE, searchWWNE_J,
                           valueWWE,  C_WWE, saveWWE,  searchWWE_J,
                           valueJE,   C_JE, saveJE,   searchJE_W,  collateralJE,
                           valueJNE,  C_JNE, saveJNE,  searchJNE_W, collateralJNE, sloanJNE, rateJNE, defaultJNE,
                           valueJEi,  C_JEi, saveJEi,  searchJEi_W, collateralJEi,
                           valueJNEi, C_JNEi, saveJNEi, searchJNEi_W, collateralJNEi, sloanJNEi, rateJNEi, defaultJNEi,
                           &capitalout, &laborout, &taxout, &taxoutSEA, genmoments,
                           distJE, distJNE, distWWNE, distWWE, distUnE, distUnNE, distUiNE, distUiE, distJEi, distJNEi);
                #endif
            #endif

            time ( &rawtime );
            time (&timeofend);
            timeinfo = localtime ( &rawtime );


            printf("===================SIMULATION ended: %s in %f seconds =================== \n",asctime(timeinfo),difftime(timeofend,timeofstart));

            //getchar();


            // Compute convergence
            critprice = fabs((laborout-lstar)/lstar);
            critprice = max(fabs((capitalout-capitalin)/capitalin), critprice);
            critprice = max(fabs((taxout-taxrate)/taxrate), critprice);

            printf("TAX: simul: %20.10f\t, old: %20.10f\t, new: %20.10f\n",taxout, taxrate, (relaxTax*taxout+(1.0-relaxTax)*taxrate));
            printf("Capital: simul: %20.10f\t, old: %20.10f\t, new: %20.10f\n",capitalout, capitalin, (relaxK*capitalout+(1.0-relaxK)*capitalin));
            printf("Labor: simul: %20.10f\t, old: %20.10f\t, new: %20.10f\n", laborout, lstar, (relaxL*laborout+(1.0-relaxL)*lstar));
            
            #if SEA == 1
                critprice=max(fabs((taxoutSEA-lumpsumSEA)/lumpsumSEA), critprice);
                printf("LumpsumSEA: simul: %20.10f\t, old: %20.10f\t, new: %20.10f\n",taxoutSEA, lumpsumSEA, (relaxTax*taxoutSEA+(1.0-relaxTax)*lumpsumSEA));
                lumpsumSEA = (relaxTax*taxoutSEA+(1.0-relaxTax)*lumpsumSEA); // if financed by laborincome tax, does not care of it.
            #endif
            
            
            // Compute Next Aggregate
            capitalin   = (relaxK*capitalout+(1.0-relaxK)*capitalin);
            lstar       = (relaxL*laborout+(1.0-relaxL)*lstar);
            
            // Update prices if allowed.
            //Rannual = alphapar*pow((lstar/capitalin),(1.0-alphapar)); // corporate sector also pay the wedgerate, so this rate is ok.
            //RquartG = pow((1.0+Rannual),0.25);    // when TFPpar = 0.25.
            //rstar   = RquartG - deltapar - 1.0 - wedgerate;             // rate r^d
            RquartG = alphapar*TFPpar*pow((lstar/capitalin),(1.0-alphapar));
            rstar   = RquartG - deltapar - wedgerate;             // rate r^d
            wstar   = (1.0-alphapar)*TFPpar*pow((capitalin/lstar),(alphapar));
            taxrate = (relaxTax*taxout+(1.0-relaxTax)*taxrate);
  
            printf("PRICES:: r=%f, rq=%f, ra=%f || wq = %f \n \n", rstar, RquartG, Rannual, wstar);
            
        } // end on BETA condition
        
        iterprice++;
#if OPT_GE == 1
    }//while(critprice > epsilonprice)
#endif



    /*_________________________________________________________________________________________________________________________________________________*/

    //////////////////////////////
    /*********** GMM ************/
    //////////////////////////////

    /****** COMPUTE THE DISTANCE ******/
    GMMmin = SMM(obsmoments, genmoments, covmat);

    /****** SAVE IT IN LOG FILE ******/
    FILE *logfile;

    logfile = fopen(tempfileout, "a");
    fprintf(logfile,"\n GMM :: %20.15f\n\n", GMMmin);
    fclose(logfile);

    #if CRS == 1
        char buffer2 [80];

        sprintf(buffer2,"INPUT/SMM_result/SMM_%f.out",GMMmin);
        strncpy(SMM_fileout, buffer2, sizeof(SMM_fileout) - 1);
        SMM_fileout[sizeof(SMM_fileout) - 1] = 0;

        logfile = fopen(SMM_fileout, "w");
        fprintf(logfile,"SMM distance = %20.15f  \n \n \n", GMMmin);

        fprintf(logfile,"------ Estimated Model parameters------\n");
        fprintf(logfile,"%20s\t\t%20.15f\n","betapar",betapar);
        fprintf(logfile,"%20s\t\t%20.15f\t%20.15f\t%20.15f\n","gmapping",gmapping[0]/meang, gmapping[1]/meang, gmapping[2]/meang);
        fprintf(logfile,"%20s\t\t%20.15f\n","KappaW",kappaW);
        fprintf(logfile,"%20s\t\t%20.15f\n","kappaE",kappaE);
        fprintf(logfile,"%20s\t\t%20.15f\n","psiparW",psiparW);
        fprintf(logfile,"%20s\t\t%20.15f\n","psiparE",psiparE);
        fprintf(logfile,"%20s\t\t%20.15f\n","nupar",nupar);
        fprintf(logfile,"%20s\t\t%20.15f\n","pz",pz);
        fprintf(logfile,"%20s\t\t%20.15f\n","sigz",sigz);
        fprintf(logfile,"%20s\t\t%20.15f\n","bankruptcy_cost",court_cost);

        fprintf(logfile,"\n");

        fprintf(logfile,"------ Estimated Moment [target] ------\n");
        fprintf(logfile,"%20s\t\t%20.15f\t\t[%20.15f]\n","KY",genmoments[0],obsmoments[0]);
        fprintf(logfile,"%20s\t\t%20.15f\t\t[%20.15f]\n","default_rate",genmoments[1],obsmoments[1]);
        fprintf(logfile,"%20s\t\t%20.15f\t\t[%20.15f]\n","E exit rate",genmoments[2],obsmoments[2]);
        fprintf(logfile,"%20s\t\t%20.15f\t\t[%20.15f]\n","E rate",genmoments[3],obsmoments[3]);
        fprintf(logfile,"%20s\t\t%20.15f\t\t[%20.15f]\n","U rate",genmoments[4],obsmoments[4]);
        fprintf(logfile,"%20s\t\t%20.15f\t\t[%20.15f]\n","U>E",genmoments[5],obsmoments[5]);
        fprintf(logfile,"%20s\t\t%20.15f\t\t[%20.15f]\n","ratio median",genmoments[6],obsmoments[6]);
        fprintf(logfile,"%20s\t\t%20.15f\t\t[%20.15f]\n","neg income",genmoments[7],obsmoments[7]);
        fprintf(logfile,"%20s\t\t%20.15f\t\t[%20.15f]\n","W>E_Q1",genmoments[8],obsmoments[8]);
        fprintf(logfile,"%20s\t\t%20.15f\t\t[%20.15f]\n","W>E_Q2",genmoments[9],obsmoments[9]);
        fprintf(logfile,"%20s\t\t%20.15f\t\t[%20.15f]\n","W>E_Q3",genmoments[10],obsmoments[10]);

        fclose(logfile);
    #endif

    //stop watch //
    time ( &rawtime );
    time (&timeofend);
    timeinfo = localtime ( &rawtime );
    tempfileoutfile=fopen(tempfileout, "a");
    setbuf ( tempfileoutfile , NULL );
    fprintf(tempfileoutfile,"\n ===================Program ended: %s in %f seconds=================\n",asctime(timeinfo),difftime(timeofend,timeofstart));
    fclose(tempfileoutfile);
    printf("\a\a\a");


    return GMMmin;


}
